Boson-Fermion pairing in Bose-Fermi mixtures on ID optical lattices 
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Boson-fermion pairing is considered in a discrete environment of bosons and fully spin-polarized 
fermions, coupled via an attractive Bose-Fermi Hubbard Hamiltonian in one dimension. The results 
of the T-matrix approximation for particles of equal mass and at double half filling are compared 
with the results of exact diagonalization and with Quantum Monte Carlo results. Satisfactory 
agreement for most quantities is found. The appearance of a stable, weak-coupling pairing mode is 
also confirmed. 
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One of the most intriguing aspects in the field of 
cold atoms involves the study of mixed Bose-Fermi 
(B-F) systems. Several boson-fermion mixtures have 
been realized, both in [l[ and without Q an optical 
lattice, and their properties have been studied. Sev- 
eral very interesting phenomena unique to mixed BF 
systems have been predicted, including for example 
the possibility of forming composite fermions through 
the pairing of a boson and a fermion 0, 0]. On 
the theoretical side, B-F mixtures have been studied 
in mean-field approximation 0, 0] and with various 
methods for treating B-F mixtures in one dimension 
(ID). This includes exact solution using the Bethe 
ansatz 0], bosonization techniques § and a Quan- 
tum Monte Carlo treatment for a B-F mixture on an 
optical lattice [?|. There remains, however, the need 
for theoretical approaches capable of reliably describ- 
ing B-F mixtures in more than one dimension [Icl |. 

A possible theoretical approach to mixed boson- 
fermion systems was proposed recently in the con- 
text of nuclear physics, to provide a framework for 
describing the transition from a fermi gas (of quarks) 
to one of composite fermions (nucleons), i.e. bound 
three-quark states. This problem was simplified by 
assuming that two of the quarks are strongly bound 
and form a boson. In this way, the extremely complex 
in-medium three-body problem was replaced by 
the much simpler two-body problem of fermion-boson 
scattering, for which a T matrix approach was devel- 
oped .12]. An interesting result of that study was that, 
due to the presence of a Fermi surface, a stable B-F 
branch was created for an arbitrarily small B-F at- 
traction. The underlying mechanism turned out to be 
analogous to the formation of stable Cooper pairs in a 
pure fermi-gas, though in the latter case the pairs are 
boson-like, whereas here they are fermion-like. Sev- 



eral interesting questions follow naturally. On the one 
hand, we would like to know how reliable the informa- 
tion provided by the T-matrix approach developed in 
[l2| is when dealing with complex systems involving 
bosons and fermions. If it is found to be acceptably 
reliable, we might then hope to further develop the 
method for application to the variety of systems in 
which boson and fermions degrees of freedom coexist, 
such as those that arise in cold atomic gases. 

As the next step in this program, we report in this 
work a study of B-F pairing in ID optical lattices 
using the above T-matrix approach. Such systems 
can be treated stastistically exactly for a fairly large 
number of bosons and fermions using Quantum Monte 
Carlo methods [H, thereby providing an appropriate 
testing ground for our method. As we will see, even 
though ID is a quite unfavorable case for the appli- 
cability of a T-matrix approach in ladder approxima- 
tion, most quantities are nevertheless reproduced rea- 
sonably well with this approach over a large range of 
coupling strengths. 

To be more specific, we consider the bosons and 
fermions on a ID lattice governed by a Hubbard model 
hamiltonian. The B-F interaction is assumed to be 
attractive and the B-B interaction to be repulsive. We 
assume further that there are no interaction among 
fermions, although a repulsive F-F interaction should 
not alter qualitatively our conclusions. The Hubbard 
model hamiltonian for such a system is given by 

<ij> <ij> 
i i 

where bi and C{ are the bosonic and fermionic annihil- 



2 



iation operators on site i, respectively. The intersite 
spacing a is set to unity, i.e. a = 1. In this work we 
carry out our analysis at zero temperature, i.e. T = 0. 
For the QMC simulations, we use a small but finite in- 
verse temperature f3 = 2L which is sufficiently close 
to the ground state. Extrapolations to larger f3 do not 
lead to significantly different results. 

In the standard T-matrix approach [13l ] to this prob- 
lem, we must solve the Bethe-Salpeter equation for the 
B-F propagator, that is [12| 

G k , k ,(K,£) = G° (K,E) <5 k , k , 

+ J2Gi(K,E)V^G klM ,(K,E) , (2) 



with G(K,E) = E k , k ' G k,k' («:,£). For a discrete 
number of sites, it is also useful to consider the equiv- 
alent diagonalization problem which can be thought 
of as a B-F Random Phase Approximation in the 
particle-particle channel in complete analogy with the 
pure fermion case [l4[. The relevant RPA eigenvalue 
equation is 
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with G k ,k' (K, -B) being the Fourier transform of the 
B-F propagator 



Here {, } denotes an anticommutator, K is the center- 
of-mass momentum of the pair and k, k' are relative 
momenta. The interactions in momentum space are 
(i/ = (&/,&&)) 



(3) 



and G k (K, E) is the B-F propagator in the HF ap- 
proximation, i.e. 
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In these expressions, Uq represents the boson conden- 
sate, p (h) refers to a fermion momentum above (be- 
low) the Fermi surface, q is a boson momentum and 

£ k ^ = — 2t^^ cos(k). The single particle energies 
in the denominator of Q are the fermion and boson 
HF energies, which are taken together with the un- 
renormalized interaction, as is usual in the HF-RPA 
self-consistent scheme [14j . Screening terms which 
would renormalize on the same footing the interaction 
and the single particle self energies are not considered 
in this prototype study but may be included in later 
work. 

Equation ([2]) can be solved analytically and the prop- 
agator summed over relative momenta is given by 
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where we define a(q) = ^/l + <5 q ,o«o an d , , = 
g v <5 P +q.p'+q'. In terms of the amplitudes X,Y the 
propagator (|3|) can be written in the following spectral 
representation: 



G kX (K,E)=J2 
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with X £+ 



(*k *£ K -k) aud xT = «f,K-k 
For a system with TV particles, the propagator has 

poles at E Pia — ±(Eq — E^ 2 ), i.e. essentially at the 
excitation energies of the A^±2 systems. The RPA am- 
plitudes obey the following normalization conditions: 



Y K Y K _ e k) + E ^P,K-p^p,K-p - ' 
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To assess the quality of our T-matrix approxima- 
tion for the in-medium B-F problem, we consider fi- 
nite numbers of sites, where exact solution is possible. 
We first restrict to half-filling and to Nb — Np and 
tb = tf. We exactly diagonalize the 6-site problem, 
i.e. Nb = Np = 3. The maximum size of the matrix 
is 318 x 318, whereas the RPA matrix has dimension 
4x4. The range of Ubf values is arbitrary. However, 
since we are interested in B-F pairing, the range of 
Ubb values is, in principle, restricted to Ubb > 1^4/1/2- 
Smaller values of Ubb lead to phase separation, i.e. the 
B-F pairs cluster together and occupy only half of the 
available space @. On the other hand, our T-matrix 
approximation does not allow [/^-values greater than 
Ubb — \Ubf | , since Ubb is only purely treated in our 
theory, i.e. only in HF approximation. Improving on 
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this point is possible but left for the future. We there- 
fore limit ourselves to \Ubf\ > Ubb > \Ubf\/2. As an 
intermediate value we take E4& = | |f/&/| throughout. 




FIG. 1: (Color online) Excitation energy as a function 
of \Ubf\ for the 6-site B-F Hubbard model and K = tt/3 
(a) and for K = (b). The solid lines correspond to the 
results of exact diagonalization and the dashed lines to the 
RPA results. 

In Fig. [TJ we present some typical examples of ex- 
cited states. We see that the agreement between ap- 
proximate and exact excitation energies is quite satis- 
factory. We should note that we have chosen examples 
where in the exact case there are only low degenera- 
cies at Ubf = 0, since RPA, because of its very low 
dimension, can not well reproduce a high degeneracy 
of uncorrelated configurations, even if there are also 
cases in RPA where in the uncorrelated limit degen- 
eracies occur. However these degeneracies are always 
less numerous than in the exact case. This is natural 
because of the dramatically reduced size of the RPA 
matrices with respect to the size of the exact ones. It 
should nevertheless be noted that the RPA excitation 
energies somehow represents the average trend of the 
bunch of exact levels. 

So far we have not invested effort to exactly diago- 
nalize problems with higher numbers of sites, since the 
dimensions of the matrices grow exponentially. How- 
ever, exact QMC results for ground state properties 
for higher number of sites are available [§] . For L = 70 
we show a comparison of the exact ground state energy 
Eq with the RPA in Fig. [2] (see 12[ for the expression 
for Eq). The error bars for QMC are smaller than the 
point size. Again the agreement is reasonable up to 
rather high values of \Ubf\- 
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FIG. 2: (Color online) Ground state energy per site as a 
function of \Ubf\ for L — 70. The solid line corresponds 
to the QMC, dotted line to the HF approximation and 
dashed line to the RPA. 



We also calculated the occupation numbers 
n b f(K) = Y^ kk ' < c K-k & k kk'CK-k' > of the B-F pairs 
and compare them with the exact results in Fig. (3^- 
In our case the B-F occupation numbers are simply 
obtained from the residue of the B-F Green's function 
© at the poles E p . We see that the agreement with 
the exact result is quite satisfactory. One might won- 
der about the upward tendency around K = ir/2 in 
the BF occupation numbers for the RPA data, since 
it is much more pronounced than in the QMC data. 
Notice, however, that for other system parameters the 
upward trend in the occupation numbers can also be 
quite pronounced even in the exact solution. Proba- 
bly, the RPA gives relatively more weight to the den- 
sity wave correlator than to the superfiuid properties 
because of the truncated model space and the fact 
that in three dimensions a gapped density wave was 
found in DMFT p]| . 

The fermion occupation numbers «/(k) =< c^Ck > 
can only be obtained in a somewhat indirect way 
within our formalism. This goes, however, completely 
parallel to what is known from the usual RPA formal- 
ism for fermions (l6| . The corresponding expression 
is given by 

n f (k) = £ |A£| 2 8( £F - 4) + E l^k P K-kl 2 e(4 - e F )(9) 

P K,p 

One can demonstrate that © conserves fermion par- 
ticle number. In Fig.[3j3 the fermion occupation num- 
bers are shown for two cases of the coupling Ubf = 
—2, —6. Again the agreement with the exact case is 
rather good, in spite of the fact that for Ubf = —6 our 
solution shows a small Fermi step whereas the exact 
one seems to be completely smooth. One should real- 
ize, however, that our theory can not describe a non- 
Fermi liquid behavior and this seems to be the case 
for stronger interactions in the exact solution. This 
deficiency of our approach to describe specific ID fea- 
tures shows up more dramatically for the boson occu- 
pation numbers where we still obtain a large fraction 
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FIG. 3: (Color online) (a)BF occupation numbers as a 
function of the total momentum K, and (b)Fermion occu- 
pation numbers as a function of the relative momentum k 
{U bb = 0.75\U bf \, L = 70). 



of particles condensed into the q = state whereas the 
exact solution is, of course, totally distributed with 
no particles in the condensate. In spite of this fail- 
ure, we conclude, however, that our first objective of 
the work has been realized, namely our T-matrix ap- 
proach which has been applied for the first time to the 
in-medium B-F case in [l2[ seems to work reasonably 
well. In addition, for t b ^ tf, our formalism is still 
valid. We compare in Table [J the ground-state ener- 
gies obtained by RPA, HF, and QMC fori/ = At b , and 
see that the results are still of the same quality as for 
the case with equal tunneling amplitude (see Fig. [J). 
Similar conclusions hold for the occupation numbers 
and for a comparison with exact diagonalization for a 
system size L = 6. We also should note that we even- 



{U bf ;U bb ) 
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(-40;6) 
(-50;10) 
(-60;12) 


-24.30 
-29.71 
-35.69 


-12.85 
-14.88 
-17.15 


-20.78 
-25.54 
-30.33 



TABLE I: Comparison of the ground state energies per 
site obtained by RPA, HF and QMC for a system of size 
L — 30 and tunneling amplitudes tf — 4t b . The HF and 
RPA relatives errors are roughly constant. 

tually would like to apply our theory to the 3D case 
where mean field and RPA theories usually perform 
much better, and where no QMC data are available. 

Another objective of this work is to confirm a sur- 
prising finding discussed in (p3 |. the appearance of 



a second stable B-F branch at arbitrary small B-F 
(attractive) interaction. It exists only because of the 
presence of a sharp Fermi surface and is thus analo- 
gous to the formation of Cooper pairs in pure Fermi 
systems. For this investigation, we pass to the con- 
tinuum limit in which the B-F propagator in HF ap- 
proximation is given analytically by 
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In Fig[4] we show the spectral function, i.e. 
—Im(G(K,E))/Tr, for a typical set of parameters for 
Nb = Njr = L/4, away from half filling. Besides the 
low lying peak which corresponds to the free fermion 
dispersion in the U b f — > limit, we see the striking 
feature that a second stable branch develops right be- 
low the continuum. This stable second branch exists 
for arbitrarily small values of U b f ■ It is the same phe- 
nomenon as was seen in our earlier work [T3|. The 
existence of a Fermi surface entails a logarithmic di- 
vergence of Re(G°) at E/2 = ep and then there is 
always a sharp state below the continuum solution 
of 1 — gbfG°(K, E) = 0. However contrary to what 
happens in the homogeneous continuum case [r2| , this 
second branch does not interact strongly with the low- 
est free fermion-like branch. In [l2j there was a level 
crossing of the two branches which does not occur 
here. Also the upper branch has rather little spec- 
tral weight compared with the lower one. It would 
be interesting to see whether this second branch can 
be found experimentally. In ID, however, we should 
note that this slightly detached second branch is cer- 
tainly an artefact of the T-matrix approximation be- 
cause the correlated fermion occupation numbers do 
not show any discontinuity (see fig. 3). However in 2D 
or 3D, we think that this second branch should exist. 
Whether it survives in a trap geometry [ItJ remains 
to be seen. 

In conclusion, we have investigated in this work 
boson-fermion pairing in a Bose-Fermi mixture on a 
ID optical lattice. The in-medium Boson-Fermion 
scattering problem was solved in T-matrix approx- 
imation. As in a previous investigation, in the 
thermodynamical limit, two stable branches were 
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FIG. 4: (Color online) (a) The solid line with circles represents the free fermion-like branch of the spectral function, the 
solid line the second B-F branch, the solid lines with squares the limits of the continuum and the dashed line a second 
plateau in the continuum, (b) The spectral representation corresponding to (JSJl for K = tt/4, see vertical broken line of 
(a) (U bf = -10, k F = 7I-/4). 



found. One corresponds to the elastic scattering 
of the fermions off the Bose condensate and the 
other is created from scattering of bosons out of the 
condensate with fermions above the Fermi sea. The 
latter comes because of the presence of a sharp Fermi 
surface and therefore has the same origin as the 
Cooper pole in a pure two component Fermi system. 
While in ID systems, the second branch is unphysical 
and just an artefact of the method, we think that in 
3D such a branch should be real. We checked the 
validity of our approach versus exact results available 
for a finite number of sites. For most quantities 
we found good qualitative and semi-quantitative 
agreement. This is satisfying because, as already 



mentioned, RPA generally works better in higher 
dimensions. More elaborate studies of B-F pairing 
are under way. 
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